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ABSTRACT 

This paper addresses the accuracy of using a 
one-dimensional analysis to predict frequencies of 
elastically-coupled highly- twisted rotor blades. Degrees 
of freedom associated with shear deformation are stat- 
ically condensed from the formulation, so the analy- 
sis uses only those degrees of freedom associated with 
classical beam theory. The effects of cross section de- 
formation (warping) are considered, and are shown to 
become significant for some types of elastic coupling. 
Improved results are demonstrated for highly- coupled 
blade structures through account of warping in a lo- 
cal cross section analysis, without explicit inclusion of 
these effects in the beam analysis. A convergence study 
is also provided which investigates the potential for im- 
proving efficiency of elastically-coupled beam analysis 
through implementation of a p- version beam finite ele- 
ment. 

INTRODUCTION 

There is a potential for improving the performance, 
aeroelastic stability, and vibration characteristics of ro- 
torcraft through the use of elastically-coupled compos- 
ite rotor blades. A complete analysis of such systems re- 
quires use of comprehensive aeroelastic codes that, due 
to their complexity and size, are in practice limited to 
modeling the elastic rotor blade with a one-dimensional 
theory. One primary function of any aeroelastic code 
is to perform an accurate dynamic analysis of the ro- 
tating blade in a vacuum. The accuracy of this part 
of the analysis may be questionable with the presence 
of elastic couplings because nonclassical beam effects, 
which are generally ignored, may become important for 
these structures. 

Many past studies have investigated the capabili- 
ties of one-dimensional analysis in predicting behav- 
ior of elastically-coupled composite beams. Smith and 
Chopra (1991), Rehfield (1990), and Nixon (1989) 
showed the importance of nonclassical effects such as 
transverse shear deformation and warping in the static 
behavior of these structures. The influences of shear 
deformation and warping in nonrotating dynamic anal- 
ysis of coupled beams have been investigated by Kos- 
matka (1988) and Kosmatka and Ie (1991). These stud- 
ies demonstrated the importance of out-of-plane shear- 
dependent warping and in-plane warping (anticlastic 
deformations) in the free-vibration analysis of beams 
in which shear deformation has significant effects, such 
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as the bending frequencies of short beams and higher 
modes of long beams. Shear deformation also is an 
important consideration for beams with some^ types 
of elastic couplings, such as bending-shear elastic cou- 
pling. Based on the work of Kosmatka and Ie (1991), 
it appears necessary to include the shear -related warp- 
ing effects for an accurate prediction of frequencies of 
bending-shear coupled beams. 

The effects of shear deformation on rotating beam 
dynamics were examined by Smith and Chopra (1992) 
in a study which extended the rotor analysis known 
as UMARC (University of Maryland Advanced Rotor 
Code, Hong and Chopra 1985a, 1985b) to include ex- 
plicit shear degrees of freedom. Results of this study 
showed improvements in the prediction of lower mode 
frequencies for bending-shear coupled beams. It is clear 
from this study that shear deformation effects must be 
included in the beam analysis to obtain accurate fre- 
quencies of bending-shear coupled beams, but the ap- 
proach of using explicit shear degrees of freedom in- 
creased the size and complexity of the formulation. 
This is undesirable in a comprehensive code because 
the additional degrees of freedom must be addressed in 
the aerodynamic models as well. 

The formulation for rotating beams is more involved 
than that for nonrotating beams because the rotation 
effects can only be included through use of the geo- 
metrically nonlinear theory of elasticity. For a general 
anisotropic beam, accounting for all the possible non- 
classical beam effects in a nonlinear formulation is un- 
desirable because of size and complexity considerations. 
One possible approach for simplifying the formulation 
is to split the equations associated with the geomet- 
rically nonlinear three-dimensional theory of elasticity 
into a nonlinear one-dimensional set of equations and a 
linear two-dimensional set of equations. This approach 
has theoretically been shown appropriate for twisted 
nonhomogeneous anisotropic blades through use of a 
variational-asymptotical method by Hodges and Atil- 
gan (1991). An explicit formulation of this approach 
is proposed which considers the influence of nonclassi- 
cal effects only on the effective beam stiffness proper- 
ties and eliminates degrees of freedom associated with 
shear deformation through static condensation. This 
formulation leads to a rotating beam analysis based on 
only those degrees of freedom which have been used 
for classical beam analyses. A second analysis, which 
should consider all the possible nonclassical influences, 
but may be based on linear theory, is used to determine 
the effective beam properties for the first analysis. If 
such an approach is accurate for geometries and materi- 
als typical of rotor blades, then rotor analyses based on 
isotropic materials and classical beam theory may be 
modified to incorporate composite materials and non- 


classical effects. 

This paper will address the accuracy of using one- 
dimensional analysis for the prediction of rotating beam 
frequencies of elastically- coup led, highly twisted rotor 
blades. There are three objectives for this study: 1) 
show that the degrees of freedom associated with shear 
deformation may be statically condensed from the anal- 
ysis, 2) show that the nonclassical influences associated 
with cross section warping, which may become signif- 
icant as a result of elastic coupling, can be accounted 
for without the incorporation of these effects explicitly 
in -the rotating beam analysis, and 3) determine the po- 
tential improvement in efficiency by using higher-order 
displacement approximations in a finite element imple- 
mentation. 

A rotating beam analysis was developed based on 
a formulation of nonlinear equations of motion and a 
finite element implementation. The formulation is de- 
rived in Appendix A to show how shear deformation 
and warping enter the theory. The formulation is non- 
linear as is required to capture the centrifugal stiffen- 
ing effects even in the linearized form of the equations. 
The degrees of freedom associated with shear defor- 
mation were eliminated through static condensation of 
the linear force-displacement relationships. The linear 
part of the formulation was implemented as a p- version 
beam finite element such that the degree of polynomial 
approximation for the bending, torsion, and axial dis- 
placements may be independently selected. This imple- 
mentation is described in Appendix B along with the 
results of a convergence study. This convergence study 
shows the efficiency of certain displacement approxima- 
tions for a bending-twist-coupled beam. 

Results of the present rotating beam analysis 
are compared with those produced by Smith and 
Chopra (1992) for a set of elastically-coupled rotor 
blades to show that static condensation of the shear de- 
grees of freedom is valid for the modes considered. At- 
tention is then focused on nonclassical effects (shear de- 
formation and warping) and their influence on the pre- 
diction of .both rotating and nonrotating frequencies for 
elastically-coupled and highly twisted beams. Compar- 
isons are made with experimental results obtained by 
Chandra (1988), 1-D analytical results obtained with 
UMARC as presented by Smith and Chopra (1992), 
and 3-D analytical results obtained using the analysis 
of Hinnant (1992). 

FORMULATION 

For this formulation, the blade is assumed to 
be a long and slender beam, and constructed from 
anisotropic materials such that displacement modes 
may be elastically coupled. The blade may deform in 
extension u t) lag bending v t flap bending w } and tor- 
sion <f > , and both built-in pretwist and elastic twist de- 
formation may be large. The equations of motion are 
formulated based on the form of Hamilton’s variational 
principal typically used in rotor analysis, 

611 = f\6U-6T-6W)dt = 0 (1) 

Jii 

The potential energy variation 8J is developed entirely 
by the elastic strain of deformation, the kinetic energy 
variation ST is developed from blade velocity terms, and 


the work variation IW is zero in the present formula- 
tion (no external loading is considered, the system is 
conservative). 

Strain Energy Formulation 

The elastic strain energy is derived in Appendix A. 
Only the linear contribution to the strain energy varia- 
tional is of interest in the present study, which is given 
by r 

RJlin = / SvikijVj dx (2) 

Jo ,f 

where (i,j = 1,9), and 

= { &i' e bo\ 6w\ &f> ty 6w‘ 6w" tv 9 6v" } (3) 

and kjj is a 9x9 cross section stiffness matrix. The 
fourth, sixth, and eighth rows and columns of k^ are 
zeros because the strain energy terms associated with 
6w* t and 6v 9 are nonlinear. The linear stiffness ma- 
trix k ij can thus be reduced to a 6x6 coupled stiffness 
matrix with diagonal stiffnesses corresponding to an ax- 
ial stiffness, two shear stiffnesses, a torsional stiffness, 
and two bending stiffnesses. For a static problem, the 
force-displacement relationship is given by: 
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where Q are forces in the directions indicated by sub- 
scripts and M are moments about directions indicated 
by subscripts. This relationship may be simplified for 
beam behavior by eliminating the shear-related degrees 
of freedom from the relationship. As was shown by 
Hodges tt al (1987) it is proper to assume the shear 
forces associated with the shear deformation are zero, 
but not the shear strains because of the presence of 
coupling terms. With Q y and Q z set to zero, the 
shear deformations may be removed through static con- 
densation. This amounts to eliminating the rows and 
columns associated with shear from the compliance ma- 
trix rather than from the stiffness matrix. The compli- 
ance matrix is formulated by inverting the 6x6 cross 
section stiffness matrix, 



and after elimination of the second and third rows and 
columns may be written as 
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The bending-related compliance terms include the flex- 
ibility associated with any shear coupling present in 
the cross section. The 4x4 compliance matrix is then 
inverted to obtain the desired 4x4 form of the fully cou- 
pled cross section stiffness matrix kj ; - , which implicitly 
includes shear deformation effects. The term k ^ is thus 
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replaced by in Eqn. 2, and the vector of continuous 
displacements is reduced to 

Svi = { 6u f e ty 6w" 6u" } ( 6 ) 

This stiffness matrix is applicable to the dynamic prob- 
lem assuming the dynamic effects associated with shear 
deformation are small. 

Kinetic Energy Formulation 

For completeness, the kinetic energy is also derived 
in Appendix A. As was the case for the strain energy, 
only the linear contribution to the kinetic energy vari- 
ational is of interest in the present study, and is given 
by 



with the vector of displacement variations for the ki- 
netic energy formulation given by, 

Stii = { 6a t 6v 6v ' £w 6w' &f> } ( 8 ) 

rriij is the mass matrix which includes rotational iner- 
tias, Cij is the linear damping matrix which is zero in 
the absence of precone (/?p)> and k,j contains the cen- 
trifugal stiffening terms which are of nonlinear origin. 

Implementation 

The linear parts of the strain and kinetic energies 
defined in Eqns. 2 and 7 were used to develop a p- 
version beam finite element so that the degree of poly- 
nomial approximation for the bending, torsion, and ax- 
ial displacements may be independently selected. Inte- 
grations over the element length were performed sym- 
bolically to increase computational efficiency of the 
analysis. Further description of the beam element for- 
mulation is provided in Appendix B. The final form of 
the rotating blade equations after application of Hamil- 
ton’s principle in discretized form is given by 

Af ij <7; + Cij qj + Kij qj = 0 (9) 

where Af i; , Cij, and Kij are the element mass, damp- 
ing, and stiffness matrices, respectively, qj represents 
the vector of discrete displacements. The elements are 
assembled to form a global system which is solved us- 
ing standard eigenvalue techniques to obtain modes and 
frequencies. 

APPLICATION 

The capabilities and limitations of the present anal- 
ysis with respect to mode and frequency predictions of 
highly-twisted el as tic ally- coupled beams are examined. 
The present analysis, referred to as CORBA (Compos- 
ite Rotating Beam Analysis) for clarity, is first verified 
for simple cases where the elastic coupling influences 
are sfnall. The predictions of CORBA are then exam- 
ined for cases where the elastic coupling effects become 
significant. Convergence of the CORBA results was 
achieved using five beam elements with cubic polynomi- 
als for the bending displacements, and quadratic poly- 
nomials for the axial and torsion displacements. These 
approximations gave convergence in the most highly 
twisted rotating beams considered in this study, and 
were more than adequate for the untwisted cases. 


Table 1: Composite blade stiffnesses for Senes 1. 


Stiffness 

Baseline 

Sym. 

Anti-Sym. 

EA/moQ J R 2 

378.1 

378.1 

378.1 

GAy/mo&R 2 

50.77 

50.43 

50.77 

GAt/mo&R? 

25.85 

25.85 

25.85 

GJ/moffR* 

.003822 

.003815 

.003796 

Ely/mo&R* 

.008345 

.008345 

.008345 

Eh/mo&R* 

.023198 

.023198 

.023198 

kn/mo&R? 

0 

-33.67 

0 

kia/mofi 2 ^ 2 

0 

- 0 

0 

kM/mofl 2 /? 3 

0 

0 

.3589 

Itn/mo&R? 

0 

0 

-.1794 

kse/mofl 2 ^ 3 

0 

0 

.1796 

L, 5 /m 0 « 2 ft< 

0 

-.001311 

0 

k46/m 0 n 2 # 

0 

0 

0 


Table 2: Frequencies for the Series 1 baseline. 


CORBA ' 

UMARC Diff." 

Pred. 

(per rev) 

(per rev) (%) 

Mode 


0.7"49 

0.747 

0.23 

1st lag 

1.147 

1.146 

0.09 

1st flap 

3.398 

3.389 

0.26 

2nd flap 

4.338 

4.315 

0.53 

2nd lag 

4.590 

4.590 

0.01 

1st tor. 

7.459 

7.416 

0.58 

3rd flap 

13.61 

13.60 

0.08 

2nd tor. 


Table 3: Frequencies for the Series 1 symmetric case. 


CORBA ‘ 
(per rev) 

UMARC 

(per rev) 

Diff. 

% 

“Fred: 

Mode 

0.749 

0.747 

0.23 

1st lag 

1.143 

1.142 

0.11 

1st flap 

3.354 

3.346 

0.25 

2nd flap 

4.338 

4.314 

0.55 

2nd lag 

4.590 

4.590 

0.01 

1st tor. 

13.63 

13.62 

0.08 

2nd tor. 


Table 4: Frequencies for the Series 1 anti-symmetric 
case. 


CORBA 

(per rev) 

UMARC 

(per rev) 

Diff. 

% 

Pred. 

Mode 

0.736 

0735 

0.08 

1st lag 

1.142 

1.141 

0.07 

1st flap 

3.344 

3.389 

1.35 

2nd flap 

4.256 

4.244 

0.29 

2nd lag 

4.367 

4.367 

0.01 

1st tor. 


Analysis Verification 

Several cases were studied to verify CORBA predic- 
tions of modes and frequencies for rotating composite 
blades. Three of the case studies are presented in this 
paper. These three configurations, referred to as Se- 
ries i, were developed by Smith and Chopra (1992) to 
investigate the effects of elastically coupled rotor blades 
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Table 5: Rotating frequencies of the Series 2 anti-symmetric case at Q = 1002 RPM. 


CORBA 

(Hz) 

UMARC 

(Hz) 

UMARC* 

(Hz) 

Experiment 

« 

CORBAt 
Diff. (%) 

UMARC*t 

Diff(%) 

Pred. 

Mode 

36.53 

36.49 

43.52 

33.6 

8.70 

29.5 

1st flap 

53.89 

53.73 

62.57 

46.6 

15.65 

34.3 

1st lag 

202.8 

202.2 

247.8 

184.0 

10.2 

34.7 

2nd flap 

336.4 

328.2 

383.6 




2nd lag 

493.6 

493.7 

493.7 




1st tor. 


t Correlation with experimental results. * ijMARC without shear deformation. 


for a soft-inplane hingeless rotor helicopter. The blade 
cross section was designed to be representative of an ac- 
tual rotor system with respect to stiffness and inertial 
properties. The main structural member of the rotor 
blade was a single cell composite box beam. The ply 
orientation of the box beam laminates was adjusted 
to produce the three configurations considered here. 
The first case is uncoupled (baseline), the second is 
extension-flap shear, flap bending- twist coupled (sym- 
metric case), and the third is bending-shear, extension- 
twist coupled (anti-symmetric case). The terms “sym- 
metric” and “anti-symmetric” refer to the orientation 
of laminates with respect to the bending axes of the 
box beam, but not to the laminates themselves. The 
individual laminates themselves are arranged in a sym- 
metric configuration for all cases. The stiffness prop- 
erties associated with each case, as reported by Smith 
and Chopra (1992), are shown in Table 1. In this ta- 
ble, EA is the axial stiffness, GA y and GA X are the 
lag and flap shear stiffnesses, GJ is the torsional stiff- 
ness, and Ely and EI Z are the flap and lag bend- 
ing stiffnesses, represents the extension-flap shear 
coupling, ki 3 the extension-lag shear coupling, kn the 
extension- twist coupling, k 25 the lag shear-flap bending 
coupling, k 36 the flap shear-lag bending coupling, k 45 
the flap bending-twist coupling, and finally k«6 the lag 
bending-twist coupling. All the stiffnesses are shown 
to be nondimensionalized by appropriate factors of mo 
the mass per unit length, fi the reference rotational 
velocity, and R the blade radius. 

The rotating natural frequencies for each case as 
predicted by two analyses, UMARC and CORBA, are 
shown in Tables 2-4. All references to “UMARC” are 
understood to mean the version which has a 19 degree- 
of- freedom shear deformable beam element, unless oth- 
erwise indicated. The difference in predictions between 
CORBA and UMARC is shown to be less than one per- 
cent for all modes except the second flap mode of the 
anti-symmetric case where the difference is 1.35 per- 
cent. Comparisons studies, not shown here, also showed 
good agreement between the two analyses for highly 
twisted blades, up to 90°. These correlations indicate 
that the present analysis has accurately captured the 
effects of rotation, twist, elastic coupling, and shear de- 
formation . 

Two more case studies, designated Series 2 , were 
examined to determine the influence of higher amounts 
of elastic coupling on the frequency predictions of 
UMARC and CORBA. The cross section geometry of 
these cases was a simple single cell box beam, without 
any nonstructural mass or secondary structure, and in 


one case the layup was arranged in an anti-symmetric 
configuration while the other was arranged in a sym- 
metric configuration. The symmetric case had a [15] 6 
layup of graphite epoxy material on the top and bot- 
tom walls while the sides had a layup of [15/-15]3. The 
anti-symmetric layup was [15]6 on top and [-15]6 on the 
bottom wall, and one side was [15]e while the other side 
was [-15]e. The box had an outside width of .953 inches 
and outside depth of .537 inches, and the specimens 
were 33.25 inches long. These cases were examined be- 
cause a set of experimental results, presented by Chan- 
dra and Chopra (1989), was available for correlation 
with the analytical predictions. 

The cross section mass and stiffness properties 
of these specimens were calculated using a two- 
dimensional analysis described in detail by Smith and 
Chopra (1990). This analysis accounts for shear defor- 
mation and the out-of-plane warping associated with 
torsion, but does not consider any other warping ef- 
fects. The mass and stiffness properties developed by 
this analysis were used as input to both UMARC and 
CORBA. 

The analytical and experimental results are listed 
in Table 5 for the anti-symmetric case and in Table 6 
for the symmetric case. The importance of including 
the shear coupling effects for the anti-symmetric case 
is demonstrated by the overly stiff predictions shown 
for UMARC* (UMARC version without shear deforma- 
tion). The frequency predictions of CORBA are shown 
to agree very well with those of UMARC in both cases. 
There is a small discrepancy in the predictions of the 
second lag modes, but this amounts to less than 4 per- 
cent. Of greater importance is the discrepancy of both 
beam analyses with respect to the experimental results. 
The correlation of CORBA with the experimental re- 
sults is shown to be poor, particularly in the lag mode, 
for both the symmetric and anti-symmetric cases. The 
error is mostly likely caused by neglecting some impor- 
tant warping terms in the cross section analysis. 

Warping Influences on the Anti- Symmetric Box Beam 

The cross section analysis employed in the verifi- 
cation studies of the last section considered only the 
out-of-plane torsion-related warping. Account of this 
warping effect gave a much more flexible and accurate 
torsional stiffness value. Analogously, the shear stiff- 
ness of the beam is also decreased by warping of the 
cross section. In this case, the majority of the effect is 
due to deformation of the cross section associated with 
shear forces both inplane (anticlastic deformation) and 
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Table 6: Rotating frequencies of the Series 2 symmetric 
case at fi = 1002 RPM. 


CORBA 

UMARC 

Exp. 

CORBAt 

Pred. 

(Hz) 

(Hz) 

(Hz) 

Diff. (%) 

Mode 

36.JRT 

307 

35.20 

4.88 

1st flap 

62.79 

62.45 

53.80 

16.7 

1st lag 

205.0 

203.0 

188.0 

9.04 

2nd flap 

392.2 

378.9 



2nd lag 

729.9 

729.2 



1st tor. 


t Correlation with experimental results. 

out-of-plane. A simplified approach for including shear- 
related warping effects in a beam is to reduce the effec- 
tive shear stiffness by a factor K which represents the 
ratio of average shear stress over the cross section to the 
shear stress at the centroid. This factor accounts for the 
near-parabolic distribution of shear stress through the 
cross section in the direction of the applied shear force, 
and is generally referred to as Timoshenko’s shear cor- 
rection factor. Since the amount of warping due to a 
shear load depends on the shape and material of the 
cross section, so does the value of K . The value of K 
was determined, using the formulas derived by Cowper 
(1966), as approximately 0.85 for the anti-symmetric 
box beam. 

The influence of the shear stiffness effect on bending 
behavior was examined for the Series 2 anti-symmetric 
box beam, but with variations of the laminate ply an- 
gles. The basic ply structure of the anti-symmetric box 
beam is [0]6 on top and one side, and [-0]e on bot- 
tom and the other side, where 0 = 15° for the baseline 
anti-symmetric configuration. The ply angle was var- 
ied from 6 = 0° to 6 = 45° for this study. The beam 
was considered non-rotating so as to isolate the elastic 
effects from the rotational effects. 

For this study, the results of CORBA were com- 
pared with those of an anisotropic 3-D p-version finite 
element analysis developed by Hinnant (1992). The 
3-D analysis used four brick elements to model the 
box beam. Convergence was achieved with ninth or- 
der polynomials for displacements along the length of 
the beam, cubic polynomials along the sides of the cross 
section walls, and linear polynomials through the thick- 
ness of each laminate. The material properties of each 
brick finite element were determined by averaging the 
material properties for each ply in the laminate over the 
laminate thickness. For cases in which the box beam 
was twisted, each brick element was twisted in a con- 
tinuous manner such that the finite element model did 
not differ from the physical model by more than one 
hundredth of an inch at any point. 

Results of the ply angle sweep for the anti- 
symmetric box beam, both with and without the shear 
correction factor applied, are illustrated in Fig. 1, 
shown as a function of error in the CORBA analysis 
with respect to the 3-D analysis. The error in the first 
bending modes is shown to increase rapidly with ply 
angle, maximizing at about 0 = 25°, and then decTease 
with ply angle. This is consistent with what might be 
expected based on the Poisson effects because the Pois- 
son’s ratio of the box beam laminates follows a similar 
trend with ply angle. The cross section warping is de- 


pendent on the Poisson’s ratio, so errors associated with 
not including all the effects of warping are expected. 
The worst error is quite significant, about 16 percent 
in the first lag mode and about 6 percent in the first 
flap mode. The error in the second and third bending 
modes is shown to be higher, with error maximizing at 
about 0 = 20° * The shear correction factor is shown 
to greatly reduce these errors, giving a very accurate 
prediction in the flap modes. 





Figure 1: Error in frequency predictions as a function 
of ply angle for the anti-symmetric box beam. 

In a second approach taken to account for all warp- 
ing influences, the lag bending stiffness was determined 
through iteration (using the CORBA analysis) as that 
required to drive the first lag bending frequency to zero 
error. The error of the second and third lag bending 
modes associated with the new lag stiffness are illus- 
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trated in Fig. 2. As shown, the error in the higher lag 
bending mo<3es is reduced, with less than five percent 
error at 6 = 30° where previously the error was in the 
10 to 25 percent range. This is an important result 
because it shows that even in cases where the warping 
effects are significant, the frequencies of higher modes 
may be accurately predicted if the same is true of the 
fundamental modes. The result is not obvious because 
the importance of direct shear effects increases at higher 
modes (beam is effectively shorter). 



Figure 2: Error in lag mode bending frequency predic- 
tions after matching fundamental mode frequencies. 

The rotating box beam with [151s layup was then 
considered with the appropriate stiffness terms as de- 
veloped in the nonrotating study. The results are shown 
in Table 7 to be greatly improved over those of Table 5, 
indicating that effects associated with rotation have a 
negligible influence on the accuracy of the frequency 
predictions. 

Table 7: Rotating frequencies of an anti-symmetric 
layup box beam at ft = 1002 RPM with refined stiffness 
properties. 


CORBA 

(Hz) 

Experiment 
- (Hz ! 

CORBA 
Error (%) 

Bred. 

Mode 

35775 

331)0 

350 

1st flap 

47.04 

46.60 

0.93 

1st lag 

190.4 

184.0 

3.46 

2nd flap 

293.4 



2nd lag 

493.6 



1st tor. 


Warping Influences on the Symmetric Box Beam 

The symmetric box beam case was also examined 
as a function of ply angle in the nonrotating configu- 
ration. For the symmetric box beam case, the shear 
is uncoupled from bending and should have little effect 
on the bending frequencies. The plots of Fig. 3 show 
that there is a dependency of the error (calculated with 
respect to the 3-D analysis results) on the ply orienta- 
tion, just as there was for the anti-symmetric case. The 
error in the prediction of the fundamental torsion mode 
(which is coupled to the flap bending mode) is shown to 
increase with ply angle to a maximum at 0 — 45° , while 
the error in the lag mode (which is decoupled from tor- 
sion and flap) maximizes at about 25°. The error in 




Figure 3: Error in frequency predictions as a function 
of ply angle for the symmetric box beam. 

the higher lag and flap modes does not follow the same 
path as the error in the fundamental lag mode with 
respect to the ply angle variations. The higher modes 
are shown to improve while the fundamental lag mode 
worsens for the ply angles above 30°. 

A new torsional stiffness was determined which gave 
a zero error in the fundamental torsion mode. The pro- 
cedure used was the same iterative procedure used pre- 
viously to obtain the improved lag stiffnesses for the 
anti-symmetric case. It was found that the improve- 
ment to the torsional stiffness drove not only the fun- 
damental torsion mode error to zero, but also drove 
the flap bending mode error to near zero because of 
the coupling between the two modes. The reverse was 
found not to be true, driving the flap bending mode 
to zero error did not correct the torsion mode error. 
Since both the fundamental torsion and flap bending 
modes could be corrected by adjusting a single stiff- 
ness value, the errors associated with the flap bending 
and torsion modes were likely from the same source, 
which was probably an alteration of the torsion-related 
warping function at the high ply angles. 

An improved lag stiffness was also calculated us- 
ing the iterative procedure. The error of the lag bend- 
ing mode is attributed to out-of-plane warping associ- 
ated with bending since this mode is decoupled from 
all other modes. . 

Application of the refined torsion and lag bending 
stiffnesses improved predictions of the higher bending 
modes as shown in Fig. 4. It is interesting that the 
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Figure 4: Error in high mode bending frequency pre- 
dictions after matching fundamental mode frequencies. 

error in the higher modes, after the corrections were 
applied, are lower at high ply angles where the beam is 
highly coupled and are worse at zero degrees ply angle 
where the beam is uncoupled. 

Effects of Large Pretwist on Nonclassical Effects 

Another important influence on composite blades is 
that of the built-in pretwist. The influence of pretwist 
could create problems for the approach of the present 
formulation because it is difficult to account for a global 
effect like pre twist in the local cross section analy- 
sis. The study of Shield (1982) illustrated the signif- 
icant influence of pretwist on cross section deforma- 
tions of bars, and the study of Kosmatka (1992) showed 
that pretwist has a significant influence on the cross 
section deformations and extension-torsion behavior of 
solid and thin-wall airfoil sections. The static behavior 
of pretwisted elastically- coupled composite beams was 
studied by Iesan (1976), Kosmatka and Dong (1991), 
and Kosmatka (1991). These studies indicate that the 
elastic-coupling and nonclassical influences of shear- 
deformation and warping can be influenced by the 
pretwist of the beam. There are no known reports to 
date, however, indicating the magnitude of the effect 
that the pretwist may have on the dynamic behavior of 
elastically-coupled beams typical of rotor blades. 

The influence of the pretwist on the nonclassical 
effects of shear deformation and warping were examined 
for the nonrotating symmetric and anti-symmetric box 
beam cases of Series S with 6 = 15°. The error of the 
CORBA predictions as compared with the 3-D results 
are shown in Fig. 5 for pretwist angles up to 90° in the 
anti-symmetric case. The change in error is small for 
the fundamental modes, with error change less than five 
percent from the untwisted case, even in the extreme 
case of 90° of pretwist. The error in the higher lag 
modes is shown to be only slightly larger, with a change 
in the error from 9 to about 16 percent in the third lag 
mode. The change in error of both the fundamental and 
higher modes, as a function of pretwist, was negligible 
for the symmetric case. 

Convergence Study 

A convergence study was performed to determine if 
use of higher order elements is beneficial when beams 
are elastically coupled. A standard h-element is defined 
for purposes of the present discussion as one with cubic 




(b) Second and third bending modes. 

Figure 5: Error in frequency predictions as a function 
of the anti-symmetric beam pre twist. 

bending shape functions and quadratic axial and tor- 
sion displacement approximations. The equivalent p- 
version element of the present formulation has p u = 1 
and p* = 1 . Since this element is common in rotor anal- 
ysis, the convergence study will consider it a baseline for 
comparison. Elements with higher order than the stan- 
dard are referenced by their addition to the displace- 
ment approximations. For example, M Std.+ lw+lt” 
refers to a beam element with one order higher approx- 
imation in flap bending and torsion than the standard 
element. 

A convergence study of a bending-twist-coupled un- 
twisted composite box beam showed slow convergence 
of the third predominantly flap mode. The cause of this 
was probably due to the coupling between the bending 
and torsion modes. Various shape function approxima- 
tion schemes were employed to determine an optimum 
for convergence of this particular mode. The results 
are illustrated in the plot of Fig. 6 which shows that 
the a Std.+lw+lt” approximation scheme had the best 
convergence. Use of that approximation scheme de- 
creased the total number of degrees of freedom from 
32 to 22, assuming a 1 percent error criteria. This 
amounts to about a one-third reduction in global de- 
grees of freedom which could relate to significant im- 
provements in run times associated with analyses of 
elastically-coupled blades. 

The composite box beam considered in the above 
study was uniform and untwisted. A second study was 
conducted on the same beam with 40° of pretwist. In 
this case, the cross section properties change as a func- 
tion of x t and, as a result, the integrations were not 
exact. Again, various shape function approximation 
schemes were employed to determine an optimum for 
convergence of the third flapwise bending mode. The 
results are illustrated in the plot of Fig. 7 which shows 
that there is no optimum. The convergence rates are 
also much shallower than those shown for the untwisted 
case in Fig. 6. This is because in addition to the elastic 
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Number of degrees of freedom 


Figure 6: Convergence of the untwisted symmetric-case 
box beam. 



The results of this study also showed that one- 
dimensional global dynamic analysis based on classical 
beam kinematics can accurately predict the bending 
and torsion frequencies of modes important to an aero- 
elastic analysis. However, the section properties used 
in the global analysis must account for the important 
nonclassical effects associated with shear deformation, 
warping, and elastic couplings. These nonclassical ef- 
fects were shown to have significant influence on the 
frequencies of the fundamental modes of highly coupled 
beam structures. Errors on the order of fifteen percent 
were reduced to less than five percent through account 
of the nonclassical effects. The influence of twist on the 
predictive capabilities of the analysis was shown to be 
small. 

The present analysis (CORBA) was implemented 
using a p- version beam finite element. Both the ad- 
vantages and disadvantages of this approach were dis- 
cussed. The p- version element proved to be convenient 
for assuring a converged solution, and allowed the de- 
sired flexibility in tailoring the displacement approxi- 
mations to the dynamic characteristics of a given beam 
configuration. Some degree of efficiency improvement 
was aemonstrated for the uniform untwisted case, but 
efficiency does not appear to be an issue for more re- 
alistic rotor blade structures. Much of the efficiency of 
using higher order elements was shown to be lost for a 
highly twisted blade. 


Number of degrees of freedom 


Figure 7: Convergence of the 40° twisted symmetric- 
case box beam. 

coupling between flap and torsion modes, the pretwist 
introduces coupling between the bending modes. The 
only higher-order element which performed well had ad- 
ditonal order increases in both bending modes as well 
as torsion. However, for this twisted case, the higher 
order elements did nothing to improve efficiency, and 
in some cases even degraded it. 

CONCLUSIONS 

A dynamic analysis has been formulated for rotating 
pretwisted composite blades which exhibit anisotropic 
behavior. The present formulation incorporated the ef- 
fects of shear deformation implicitly through elimina- 
tion of the shear variables in the material compliance 
matrix. Results showed that this approach was able 
to capture the most significant effect of shear defor- 
mation, namely the reduction in effective bending stiff- 
ness that occurs when a substantial amount of bending- 
shear coupling is present in a beam. The difference 
between implicit and explicit use of shear degrees of 
freedom was shown to be less than 2 percent up to 
the second bending modes of some representative rotor 
blades, and less than 4 percent up to the second bend- 
ing modes of some highly coupled box beam specimens. 
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APPENDIX A 

STRAIN AND KINETIC ENERGY 
FORMULATION 

As the formulation presented here is nonlinear and 
explicit, the number of terms in the energy expressions 
can quickly grow to an unmanageable size. Further, 
many of the terms may be negligible compared to other 
important terms. To reduce the number of terms to 
only those of significance, an ordering scheme is em- 
ployed where terms of 0(e n+2 ) and higher are elimi- 
nated in the presence of terms of 0(c n ). All displace- 
ment variables defined in this formulation are assigned 
an order of c with two exceptions. The axial displace- 
ment it c is of order e 2 and the twist deformation <t> is 
of order one. The latter exception results from making 
the analysis accurate for rotor blades with very large 
elastic couplings associated with twist deformation. 


Geometry and Coordinates 


The present formulation requires six coordinate sys- 
tems. Shaft , gimbalhd , and hub-fixed coordinate sys- 
tems are defined as illustrated in Fig. 8. For purposes of 



Figure 8: Geometry of the shaft, gimbal and hub. 

the present formulation, the shaft reference frame is an 
inertial system. A blade reference frame (e r ,e y ,e z ) is de- 
fined with its x-axis directed along the elastic axis of the 
undeformed blade as shown in Fig. 9. The elastic part 
of the blade is offset from the center of rotation a dis- 
tance h x i T . A cross-section reference frame (e Xi e v ,i{) 
is defined with origin at some position (h x + x)e x along 
the elastic axis of the blade, and with origin on that 
axis acting as the reference point for the cross section. 
The unit vector e v is directed along the chord direc- 
tion of the blade cross section while is defined by 
the cross product of e x and e^. Thus, the cross-section 
system is an orthonormal vector set which rotates with 
the built-in pretwist of the undeformed blade. A de- 
formed reference frame is identical to the 

cross-section set before deformation, but translates and 
rotates with the bending and twist of the rigid cross 
section plane to a new position after deformation. 

The unit vector triads of each coordinate system are 


related by the following equations: 



is 

Js 

k s 

Ig 

Jg 

k G 


i 

} 


Ih 

Jh 

Kh 



( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 


and the transformations matrices are themselves given 
by: 


[Tas] 


Phg] 


[Tbh] 


[ t cb] 


■ cos Pac o —sinflac 
0 1 0 
sin fee 0 cos 0gc 

’10 0 
0 cos 0qs sin Pgs 
0 —sin Pas cos 0gs 

cos ip sinV 1 0 
— sin V cos tl> 0 
0 0 1 . 

cos Pp 0 —sin(3p 
0 1 0 
sin Pp 0 cos Pp 

’ 1 0 0 1 
0 cos f3 sin /9 
0 — sin/9 cos/9 


(15) 

(16) 

(17) 

(18) 


The transformation between the deformed and cross- 
section systems [Toe] is derived later in this appendix. 

Strain Energy Derivation 


Consider the position of a point on the cross section 
of a rotor blade before deformation with position vector 
given by 

To = (hx + X)e x + T)e„ + (e ( (19) 

After deformation, the position vector is given by 

R = Rq -f Re + Rw (20) 

where Rq represents deformed position of the cross sec- 
tion reference point, Re represents deformation asso- 
ciated with the rigid rotation of the cross section, and 
Rw represents deformation associated with warping of 
the cross section. The position vectors are defined as 
follows: 

Ro 
Re 
Rw 


= (fc r + x + u 0 )e x + v 0 e^ + u> 0 cc (21) 
= 0£* + »?£„+<:£( (22) 
= + + M c (23) 
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Figure 9: Geometry of the elastic blade. 


t 


where W Uf and Wq are warping displacements de- 
fined as 


W u 

= «Vu/l + «>C. Q ( + + 



V 'c,^Qn + 0Cl/'uAf < + ^VuT 

(24) 


= uViM + v'c.'PnQi + O’l'I’fiMi + 



+ e C^nM( + ^VijT 

(25) 

w c 

= + 6 v ll> CM, + 



^e.V’C Q, + 

(26) 


where subscript s denotes shear strains due to shear 
deformation and 0 is rotation due to bending. The 
warping terms represent nonclassical contributions to 
the displacements as a result of cross section deforma- 
tion. The notation for the warping tf>ij gives the dis- 
placement in the direction i associated with a load j } 
and the magnitude of the displacement in the i direc- 
tion is shown to be proportional to the displacement 
associated with the load direction. The displacements 
associated with warping are in general small for beam 
structures, with only a few exceptions. The most well- 
known exception is the out-of-plane warping associated 
with torsion of noncircular beams in the present 
formulation). With a completely general approach to 
anisotropic beam theory, any of the 18 warping terms 
shown above could be significant for a particular con- 
figuration. Thus, for the general approach, all of the 
warping terms would be maintained within the order- 
ing scheme, even though for most practical cases all but 
a few terms could be eliminated. 

However, the warping displacements have little use 
in a one-dimensional analysis because the primary ob- 
jective is accurate assessment of global behavior. The 
important contribution of warping has been shown in 
past studies to be a reduction in the effective beam stiff- 
nesses. As such, the warping unnecessarily complicates 
development of the one-dimensional analysis and will 
be eliminated except for some key terms which have 
been shown to be important, even for isotropic beams. 
The other effects of warping can be captured in a de- 
tailed cross section (local) analysis which is uncoupled 
from the beam (global) analysis. 

The warping terms which are retained are the out- 
of-plane torsion-related warping V’uT , and the two out- 
of-plane shear-related warping terms and *PuQ n - 

If the Timoshenko-type shear deformation model is ap- 


plied fthe cross section is assumed to remain plane), 
then tp U Q < = C and = *7* The deformed position 

vector is then rewritten with R^, = [<£VuT + + 

<.C]£* as 


R 


({h x + x + u 0 , w 0 , u>o} + {(^VuT + 


v c.v + K.O^XWdc]) 



(27) 


where Tpc is the transformation matrix between the 
deformed and cross-section coordinate systems, and 
will be derived in the next paragraph. 

The sequence of rotations for transformation from 
the undeformed cross-section axis system to the de- 
formed axis system is {6(, -0,, where 6$ is the Eu- 
ler bending rotation in the lead-lag plane (given no 
pretwist), $ v is the Euler bending rotation in the flap- 
wise plane (given no pretwist), and 4> is the elastic twist 
which may be a large angle. The transformation matrix 
is then defined as 


Pi>c] 



‘ 1 

0 

o ■ 

■ 1 

0 

-Or, I 


0 

cos <f> 

sin <f> 

0 

1 

0 


0 

— sin <6 

COS <f> 

K 

0 

1 


i e, o i 
-0 ( 1 o 
o o 1 


(28) 


where the small angle assumption has been employed 
for the bending rotations. The rotations may be written 
in terms of the cross-section kinematic variables as 


= K,r - (29) 

= {u>c,x +Art> c )c ( (30) 

which, when substituted into Eqn. 28, gives the trans- 
formation matrix as 


[Tdc] = 


1 

v'c - 

K + vjT 

-(V' c -W e /}')C06 4> 

Pc sin^ 

sin <f> 

~(vr c + v e 0') sin <j> 

+ C0S<£ 


-(«£ + v e /?')cos^ 

Pc COS <f> 

cos <j> 

.+(«' -w e 0')sm<f> 

— 8\n<f> 

m 
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where, 


Pc = —(v'c — w e p')(w , e + v e f?) (32) 

This transformation agrees with that of Kosmatka 
(1986) if 4 is assumed to be a small angle. 

The strains are developed in terms of the displace- 
ments by substituting the derivatives of the position 
vectors into the strain component definitions as given 
in Wempner (1981). The position vectors have been 
defined in terms of the cross-section coordinates, and 
the derivatives were calculated as follows: 


r, x 

= {1,0,0} 

(33) 

r.n 

= {0,1,0} 

(34) 

f X 

= {0,0,1} 

(35) 

Rx 

= {G*i, G*2 ,G*s} 

(36) 

R,n 

= {G*)l,G ! «(2,G ; i)3} 

(37) 

R< 

= {G(1,G(2,G(3} 

(38) 


where the G, terms are defined within the ordering 
scheme as: 

G x i = 

1 + «'c - “ C*c + 

4* KCt/ c + Wc - V v cP f - (w c p f ) cos <t> 


+ (n v c + Cv'c + C*>cP' - rjw c /3‘) sin 4] 

+(*V«t)' (39) 

G x 2 = 

v' — Weft* — 4'[( cos 4> 4* rjsin 4] (40) 

G x 3 = 

w f c + V *P> + 4 , [v cos 4 - C sin (41) 

Gjji = 

“ ( v e ~ w cfi) cos 4- (w' + v e P*) sin 4 + 

(42) 

G„ 2 = 

cos 4 4“ {VcWcP* 2 - Vc^'p* + 

WcW^p 1 — u;' ) sin <£ (43) 

G»j3 = 

sin 4 (44) 

G C1 = 


and is the curvature in the flapwise plane and 
is curvature in the lead-lag plane. The strain compo- 
nent definitions simplify after substitution of the unde- 
formed position vectors to 



= (Rr -Rx- 1)/2 

(50) 

*xtl 

= (R*R,v) 

(51) 


= (R,*R,<) 

(52) 

e »7*7 

W K«I|(«0 

(53) 


where € rrj and c X ( are the engineering form of the shear 
strains. The three nonzero strains are calculated by 
carrying out the dot products. In terms of the dis- 
placements defined in the cross-section system, these 
are shown after application of the ordering scheme. 

txx = u'c + ^(v'c - W e ^) J + i(«/' + tie/?') 2 - 

V K n ~ C*c + |(»? 2 + CV 2 + (^V«t)'(5 4) 

= < + (V’uT.i) - 0<t>' (55) 

€ x< = w e. + (V’uT.C + V)4>' (56) 

These strains are defined in terms of the blade coordi- 
nate system through use of the transformation [Tcb] 
as 



= v'+l^ + l ti^ + itf+cV 2 



-t;"[»7Cos(/? + ^) - C 8 >n(/? + ^)] 



— ui"[^sin(^ -f (/>) + C sin(/? + </>)] 



+(^'V’ut)' 

(57) 

€ XfJ 

= w' cos(/? + ^) + tn' sin(/? + ^) + 



(tfruT.i) - C)4>' 

(58) 


= tv', cob(P + <t>) - v', sin(^ + </>) + 



(*i>uTX + W 

(59) 


At this point a variable substitution is made which 
eliminates the kinematic contribution of forshortening 
from the axial displacement. It has been shown by 
Kaza and Kvaternik (1977) that this substitution pro- 
vides the convenience of developing centrifugal stiffen- 
ing terms associated with forshortening in the kinetic 
energy formulation rather than the strain energy for- 
mulation. The substitution is 


w f Ct - (u/' + v c p f ) cos 4 + (v f c - W c p f ) sin 4 + 

4 f 4uT£ (45) 

G( 2 = 

- sin 4 + ( v c w c p n - v c v' c P ' + 

w c w' c ff - v* c w f c ) cos 4 (46) 

G( 3 = 

cos 4 (47) 

and the curvatures are given by 

K, = («e - I v c P" - 2w' e P' - Vj3 ,2 ) cos <j> + 

(w" + v c p" + 2v' c (3' - w e /?' 2 ) sin <f> (48) 

K C = ( w e + v cP" + 2t/'/?' ~ WcP 12 ) COS <f> - 

(v" - w e p" - 2w' e p l - v e /?' 2 )sin </> (49) 


ft/ 2 ( 60 ) 

where u' represents the elastic axial strain without 
kinematic contributions from transverse bending dis- 
placements. The strain components then become in 
final form: 

<« = ti' + ^ (t] 3 + <V 2 + (^Vut)' 

-v"[t) cos(/3 +</>)- c sin(Z? + </>)] 
—w"[T)sin(p + <t>) + C sin(/? + <f)] (61) 

c«j = cos(/? + <f>) + tv', sin(/? + 4) + 

WuT.n-O# (62) 

C X( = ti>: cob(P + 4 >)-v', sin(/? + <}>) + 

(V\.7\< + V)<t>' (63) 
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The elastic stress-strain relationships employed in this fty = 0. After application of the ordering scheme, the 
formulation are given as velocity is given by 


{£H 


8 * fr 


Q 15 

QU 


^55 

qI 





(72) 


where the Q- represent the material stiffness at a loc- 
tion in the cross section. The material stiffnesses are an 
average value based on the individual ply material and 
orientation, and also depend on the orientation of the 
laminate with respect to the cross section axes. The 
variation of the elastic strain energy can then be writ- 
ten as 


SJ = a ^xx "1" &xr)fexrf "I" ^x( ^x( } df)d(^dx (65) 

The stress-strain relations are substituted into the 
strain energy variational, followed by a second substitu- 
tion of the strain-displacement relations (Eqns. 61-63) 
for the strains. After integrating over the area, the 
strain energy variation becomes 


SU = f dviVijij dx + f SbiDidx 
Jo Jo 


( 66 ) 


where (i } j = 1,9), and 


K = 

u + u>ft y — vft, — (ft, -f v f )rj cos 02 + 

(fty - w'X cos + *?fty sin 02 + ft*C sin 02 + 

C COB 02V* -I- T) sin 02 V*4> — T)CO8 02tV f 4> — 

C sin 02W f j> + 7 } sin -|- £ sin 0 7 v f (73) 

V v = 

V — tuft r -f (A* + x + ti)ft, — ft, 7}V f cos 02 — 
SlxCtOB0 2 — COS 02 ~ ft,(tE/' COS 02 ~ 

ft*rj sin 02 — rjj> sin 02 — ft, rj w ' sin 02 + 
ft,Cv'sin^2 (74) 

V, = 

ti; -|- vftx — ( h x + i + u)ft y + ft r 17 cos 02 + 

fty T)V f COS 02 + T)j> COS 02 + fty£t/ COS 02 + 

fty 7 w* sin ^2 - sin 02 — ftxC sin & — 

fty<t/sin/? 2 (75) 


Svi = { 6u' e 6v' $ Sw f s 6j> W tut 6w" bo 1 Su " } (67) 

The first integral of Eqn. 66 represents the linear part 
of the strain energy and the second term represents the 
nonlinear contribution to the strain energy. 

Kinetic Energy Derivation 

Now, the kinetic energy is derived. The position of 
a point on the deformed blade as given by Eqn. 27 may 
be written using the blade reference displacements and 
neglecting the warping displacements as 


R = ({>»* + * + U, V, w} + {o, T), C }pbB]) 



( 68 ) 


where \Tdb] is the transformation between the de- 
formed and blade coordinate systems which is given 
by 

[ Tdb ] = [Tdc][Tcb] (69) 


The velocity of a point on the deformed blade is written 
as 

a n 

V=^ + nxR (70) 

where 


a 





{O.O.ftoMTHGfPflff] 7, 



(71) 


and ft 0 is the rotation rate at which the hub spins about 
the gimballed zq axis. If there is no precone then ft r = 


where 02 = 0 + <f>- After taking the variation of the 
velocity, the following substitutions, which are based 
on Eqn. 60, are made into V and 6V. 

ti = ii € — f ( v'v ' + tt/ti/) dx (76) 

Jo 

6a = 6u e — f ( v’6v ' - w'bw') dx ( 77 ) 

Jo 

The variation of the blade kinetic energy is given by 

6T = Hi pV- SVdrjdCdx (78) 

where p is the mass density of the blade. After sub- 
stituting the velocity as defined in Eqn. 72 into the 
kinetic energy expression, calculating the velocity vari- 
ation, and carrying out the dot product, the variation 
of the kinetic energy may be written as 

R 

m{[T u ], 6ui + [Xu]i biii + T a & 2 + Tp } dx (79) 

where (i = 1 , 6 ) and the vector of displacement varia- 
tions for the kinetic energy formulation is given by, 

6u{ = { 6u e 6v 6v* 6w 6a/ 6$ } (80) 

The quantities [Ty]* and [Tyh represent groups of terms 
which may be functions of both ti and ti. The terms 
T a and Tp represent additional terms in the kinetic 
energy which result from the variable substitutions de- 
fined in Eqns. 76 and 77. T a represents the nonlinear 
anti-symmetric contribution of Coriolis force and is de- 
fined as 

T a = 2 [ (t/t/ + ii/ti;'K (81) 

Jo 
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The term 7> gives the contribution from the integral 
part of the variation substitution for 6u ey 


Tp — — (x + 2v) f (v'Sv* -h w f 6w f ) d£ 

Jo 


(82) 


After integrating by parts, 7> may be written in a more 
convenient manner as 


T f = -(F a + F cor )( v’6v' + w'6u>') 


where 


Fa 


-j: 

■f 


mx d£ 


2 mi) d£ 


(83) 

(84) 

(85) 


The terms associated with F A reflect the centrifugal 
stiffening effects on the flap and lag equations while the 
terms associated with F cor reflect the nonlinear Coriolis 
damping effects in those equations. The terms associ- 
ated with F a and F cor are added to T v • and T w # which 
allows the linear contribution to the kinetic energy vari- 
ation to be written as 

r* . 

ST Iin = / {SCiimijUj + foiCijUj + faikijUj) dx (86) 

Jo 

A more useful form of the above expression is obtained 
by integrating the variation in kinetic energy by parts 
over time. This can be done because in applying Hamil- 
ton’s principle the variation in kinetic energy will be in- 
tegrated in time. By temporarily switching the order of 
integration, the integration by parts can be performed. 

I I SuifTiijUjdx dt — I I SuitnijUjdt dx = 

Jt\ Jo Jo Jt | 


APPENDIX B 
IMPLEMENTATION 

Introduction 

The present formulation is implemented as a beam fi- 
nite element. Many past analyses for rotating blades 
have used this approach, but the order of polynomi- 
als used to approximate the displacements has varied. 
The analysis of Kosmatka (1986) uses a quadratic tor- 
sion and axial approximation along with cubic Hermi- 
tian polynomials for bending. This set of assumptions 
provides the same level of accuracy in the torsion and 
axial deformations as in the bending deformations. The 
analyses of Hong and Chopra (1985a, 1985b) and Smith 
and Chopra (1991) use similar displacement polynomi- 
als, but with a cubic axial approximation, developed as 
a mean for improving the axial mode predictions. 

A higher order element capability was developed for 
the dynamic analysis of beams in the GRASP code 
(Hodges et. a/., 1990). In this code the user could 
independently increase the order of polynomial ap- 
proximation of each displacement to match the physi- 
cal characteristics of the beam. This is the so-called 
p- version finite element approach, and seems ideally 
suited for application to analysis of elastically-coupled 
beams because of the dramatic influence elastic cou- 
plings have on beam flexibility in some displacement 
modes. The study of Hinnant (1989) demonstrated 
that, given proper modeling of the beam geometry, 
there is also substantial savings to be gained by use of 
p- version elements in terms of total number of degrees 
of freedom required to obtain an accurate solution. 

Finite Element Implementation 

The linear parts of the strain and kinetic energies 
as defined in Eqns. 66 and 88 are used to develop a p- 
version beam finite element. The continuous displace- 
ments which appear in these expressions are u, v, w t 
and and are functions of both x and time. The con- 
tinuous problem is discretized by introducing discrete 
degrees of freedom $ which are related to the continu- 
ous displacements according to 


JO Jt 1 

— f f dz dt (87) 

Jt x Jo 

U — 

t=i 

p* 

(89) 

After a similar operation on the damping term of 
Eqn. 86, the linear variation of kinetic energy becomes 

V = 

* = 1 

(90) 

f R 

6Tu„ = / fui{mijUj + CijUj + kijUj) dx (88) 

Jo 

w = 

*=i 

(91) 


4> = 

i = 1 

(92) 


where N% are shape functions defined later in 
this section. Substitution of these equations into 
Eqns. 66 and 88 gives the strain and kinetic energies 
in terms of the discrete degrees of freedom. The virtual 
energy expression defined in Eqn. 1 may also be written 
in discretized form as 


SU 


-i: 


N 


- ST,) 


«=i 


dt 


(93) 
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where N is the number of spatial elements used to dis- where = N ’f = Nf . The shape functions N? and 
cretize the elastic blade. Each element is represented jq™ require C 1 -type continuity. These shape functions 
using the discrete displacements as are given by: 


SU - ST — SqJ {Mijqj + C ijqj + K ijqj } (94) 

where the element mass, damping, and stiffness matri- 
ces are defined by 


Mij = 

[ BitmkiBijdx 

Jo 

f R 

(95) 

c fi = 

/ BikCk\Bij dx 
Jo 

f R 

(96) 

Kij = 

/ (Aik^ktAtj dx - BikktiBtj)dx 
Jo 

(97) 


where B ik = Bj- and Aik = Ajj. B is a matrix of shape 
functions and shape function derivatives which satisfies 
the relationship 

Ui = (Dr)ij[uj] = (DT)ij[Hj k qk] = Bikqk (98) 

where Uj is a vector of the continuous degrees of free- 
dom u, v, u;, and <j>. Dt is a matrix of derivative opera- 
tors associated with the kinetic energy formulation and 
H is a matrix of shape functions whose arrangement 
depends on the selection of discrete variables in g, and 
satisfies Eqns. 89-92. The definition of A ^ is similar to 
that of Bij except that it is associated with the strain 
energy formulation. Thus, B may be replaced by A and 
subscripts of T may be replaced by V in Eqn. 98. 

The discrete degrees of freedom are divided into 
two sets, external and internal. There are twelve ex- 
ternal degrees of freedom which have physical signif- 
icance as the displacements and rotations associated 
with the ends of the beam finite element (six on each 
end). These deformations are depicted in Fig. 10. The 
shape functions for N- 1 and Nf are identical and have 
C°-type continuity. There are two well-known linear 



Figure 10: Beam element showing external discrete de- 
grees of freedom. 

polynomials used to define this set: 

K = l-y (99) 

N° = y (100) 



X X 

= 2 7T- 3 7T + 1 

(101) 


x 3 ,‘ 2 

~ p 2 / + 1 

(102) 

*3 

X 3 X 2 

= -273+375- 

(103) 

Nl 

X 3 X 2 

~ / 2 l 

(104) 


where N* = N™ = N>. 

The internal degrees of freedom have no physical 
significance, but are simply coefficients of the higher 
order shape functions. The internal degrees of freedom 
serve to increase the accuracy of the transformation 
from the discrete problem having a finite number of de- 
grees of freedom to the continuous problem having an 
infinite number of degrees of freedom. In the present 
formulation, the number of internal degrees of freedom 
is limited to four for the C°-type displacements, and to 
two for the C l -type displacements. There are, there- 
fore, a total of six internal shape functions associated 
with each continuous displacement u, v, w, and <f>. The 
additional C°-type shape functions for u and <f> are 


N° = VZ(£- X j) 

(105) 

K = A(-2^ + 3^-^) 

(106) 

= ^(5y?-io£ + e£-i) 

(107) 

*n ^£ 5 „ „ X 4 ^ X 3 


*6° = -42^ + 105-^-90^- + 


_ a £ 2 n X 

30p-3 7 

(108) 


These shape functions are derived by Hinnant (1989) 
based on satisfaction of two requirements: first, the 
higher order shape functions must be zero at the ele- 
ment boundaries, and second, they must be orthogonal 
with respect to their first derivative. The additional 
C^-type shape functions for v and xv are given by 

"l = + ( 109 ) 

"l « ' / i<-7r+ 5 fj7- 2 F + l* ) <110) 

The derivation of these higher-order polynomials is sim- 
ilar to that of the C^-type polynomials, only the func- 
tions must also have zero slope at the element bound- 
aries, and must be orthogonal in their second deriva- 
tive. 

The arrangement of shape functions in the matrix 
of shape functions H depends on the arrangement of 
discrete degrees of freedom in q. To facilitate the el- 
ement assembly process, the discrete unknowns were 
grouped with the first twelve external nodes together, 
followed by the twelve internal nodes (4u, 2v, 2u;, and 
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4^). The arrangement of the vector of discrete degrees 
of freedom is given as 


* = < 


w x 

<h 


Yi 


v 


i 
«2 
V 2 
W 2 

<t> 2 

— u/ 7 

v 2 

^4 

«6 
V3 
V 4 
W3 
W 4 
<t> 3 
<^4 
<*5 


( 111 ) 


Before the symbolic integrations of Eqns. 95-97 can 
be carried out, the mass, damping, and stiffness cross 
section matrices (m, c, k> and k) must be defined as 
polynomials in x. The cross section terms are func- 
tions of x because of the presence of the twist angle in 
many of the terms, which is itself a function of x. In the 
present formulation, it is desired to have the capabil- 
ity of accounting for changes in cross section properties 
beyond that due to twist, such as taper, for example, 
A beam element does not allow for such effects directly, 
so a quadratic polynomial curve fit was adapted to in- 
crease the accuracy of the element for changes in cross 
section properties along its length. 

The mass, damping, and stiffness matrices as given 
by Eqns. 95-97 were symbolically integrated to obtain 
24 x 24 element matrices. These matrices were im- 
plemented in an analysis to determine the modes and 
frequencies of highly-twisted elastically-coupled rotor 
blades. As part of this implementation, the displace- 
ment approximations could be chosen for each contin- 
uous displacement independently. The external dis- 
placements represent the minimum number of degrees 
of freedom for each element, while the maximum is 
given by use of all twelve internal degrees of freedom. 
Any choice between 12 and 24 degrees of freedom per 
element could be accomodated in the analysis. The 
notation adopted for the present formulation is to se- 
lect a “p” value which represents the number of in- 
ternal degrees of freedom associated with a particular 
displacement. For example, an element with p u = 1 
and pj, = 1 uses the basic cubic hermitian polynomial 
approximation in bending (no internal degrees of free- 
dom) and quadratic polynomial approximations in the 
axial and torsion displacements. This particular exam- 
ple happens to represent the mo6t common approxima- 
tion used in finite element rotor blade dynamic analysis 
because it gives an equivalent level of approximation in 
all displacement modes. 
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